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ABSTRACT 

With the ever increasing size and complexity of fully self-consistent simulations of galaxy 
formation within the framework of the cosmic web, the demands upon object finders for these 
simulations has simultaneously grown. To this extent we initiated the Halo Finder Comparison 
Project that gathered together all the ex perts in the field an d has so far led to two comparison 
papers, one for dark matter field haloes ( IKnebe et al and one for dark matter subhaloes 

dOnions et al.ll2012h . However, as state-of-the-art simulation codes are perfectly capable of 
not only following the formation and evolution of dark matter but also account for baryonic 
physics, i.e. gas hydrodynamics, star formation, stellar feedback, etc., object finders should 
also be capable of taking these additional physical processes into consideration. Here we 
report - for the first time - on a comparison of codes as applied to the Constrained Local 
UniversE Simulation (CLUES) of the formation of the Local Group which incorporates much 
of the physics relevant for galaxy formation. We compare both the properties of the three 
main galaxies in the simulation (representing the Milky Way, Andromeda, and M33) as well 
as their satellite populations for a variety of halo finders ranging from phase-space to velocity- 
space to spherical overdensity based codes, including also a mere baryonic object finder We 
obtain agreement amongst codes comparable to (if not better than) our previous comparisons 
- at least for the total, dark, and stellar components of the objects. However, the diffuse gas 
content of the haloes shows great disparity, especially for low-mass satellite galaxies. This 
is primarily due to differences in the treatment of the thermal energy during the unbinding 
procedure. We acknowledge that the handling of gas in halo finders is something that needs to 
be dealt with carefully, and the precise treatment may depend sensitively upon the scientific 
problem being studied. 

Key words: methods: TV-body simulations - galaxies: haloes - galaxies: evolution - cosmol- 
ogy: theory - dark matter 



1 INTRODUCTION 

In a series of precursory papers jKnebe et al.ll201 fl : lOnions et al.l 
I2OI2I see also Knebe et al., in prep.) we have compared and 
quantified the differences arising from the application of various 
halo finders to the same (dark matter only) data sets. The moti- 
vation for such comparisons is obvious: while past decades have 
seen great strides to better understand disparities in cosmologi- 
* E-mail: alexander.knebe@uam.es cal simulations stemming from different gravity (and hydrodynam- 
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ics) s olvers (e.g.'prenk et al.'l999l;lKnebe et al.'2000';'O'Shea et ajj 
2005: Agertz et al. 2007; Heitma nn et al. 2008; Tasker et al. 20081; 
Robertson et al.l l20ld ; ISpringell Hoiol ; IVazza et al.l |2011|) it was 
about time to initiate a similar project for object finders. To put 
it in a provocative way, if you are not only interested in the matter 
distribution in the universe, any simulation is only as good as the 
halo catalogue derived from it, especially at a time when we will 
have access to extremely large redshift surveys (just to name a few, 
BOSS, WiggleZ, eBOSS, BigBOSS, DESpec, PanSTARRS, DES, 
HSC, Euclid, WFIRST, etc.) against which the next generation of 
cosmological simulations (and the haloes/galaxies found in them) 
should and will be compared against. 

However, all our previous comparisons solely focused on dark 
matter haloes and left aside any possible influence of the baryons. 
In this work we extend the comparison project to a cosmological 
simulation that not only models gravity, but simultaneously fol- 
lows the evolution of the baryonic material by incorporating a self- 
consistent solution to the hydrodynamics. The simulation further 
includes star formation and stellar feedback, all explained in more 
detail in the subsequent Section[2] Please note that we are not try- 
ing to address the question of (disk) galaxy formation in this work, 
but rather aim at using one of these simulations to verify whether a 
halo finder applied to it will find the same (galactic) objects as any 
other halo finder. In particular, we will study both the larger (disk) 
galaxies as well as their respective satellite populations. To this ex- 
tent we use a constrained zoom simulation of the Local Group of 
galaxies consisting of objects akin to the Milky Way (MW), the 
Andromeda galaxy (M3I), and M33. This data forms part of the 
Constrained Local Universe (CLUES) projectQ 

The questions to be addressed are plain and simple: 
while it is well established that baryonic physics alters the 
particulars of dark ma t ter haloes and subhalo populations 
telumenthal et al.| I986I; Tissera & Dominguez-Tenreird 1 9981; 
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2OI2I ; iDi Cintio et alj|2012h . there remains the question of how 
this will influence the performance of an object finder. Note that 
cold baryonic material tends to cluster more strongly than its dark 
matter counterpart due to the dissipative nature of the cooling 
process allowing for the disposal of energy via radiation. While star 
particles - generated on-the-fly during the simulation according to 
a physically motivated prescription that transforms cold gas into 
stars - only feel gravitational forces (like the dark matter particles), 
they nevertheless originate from compact, condensed regions. 
Additionally the gas particles themselves carry not only kinetic but 
also thermal energy giving rise to thermal pressure in the medium 
(specified by the adopted equation of state). At present, halo 
finders deal with these subtleties differently, with some properly 
accounting for the gas pressure and others ignoring it. Therefore, 
the fundamental question is, what will different halo/galaxy finders 
find when it comes to simulations including baryonic physics? 

Here we compare a range of different (halo) finding tech- 
niques ranging from phase-space to velocity-space to spherical 
overdensity based codes, but also including a mere baryonic object 
finder. This naturally raises the question about the right approach 
to actually find galaxies. Should one hunt for dark matter haloes 



http : / /www . clues-pro ject . org 



first and then characterize the galaxy residing in it? Or is it more 
reasonable to initially locate galaxies and then identify their sur- 
rounding dark matter haloes? Our suite of object finders contains 
both methods of operation and we will hence learn about these dif- 
ferent strategies, too. 



2 THE DATA 

We us e th e same simulati on already p resented in iLibeskind et al] 



ll201Qf); iGotflober et al.l feoioh; ILib eskind et al.' ('201l|); 



Kne be etal. ( 2010); Libeskind et aU (I20IL1; Knebe et al. (20lJ); 
di C intio et alj I2OI1I) ; iDaval & Libeskindl j20I2l) ; iDi Cintio et al] 
1 120 121) and refer the reader to these papers for a more exhaustive 
discussion and presentation of our constrained simulations of the 
Local Group that form part of the aforementioned CLUES project. 
We briefly summarise their main properties here for clarity. 



2.1 Constrained Simulations of the Local Group 

We choose to run our simulations using standard ACDM ini- 
tial c onditions, that assume a WMAP3 cosmology jSpergel et al.l 
|2007|) . i.e. n„ = 0.24, Sit = 0.042, Q.k = 0.76. We use a 
normalisation of erg = 0.73 and an = 0.95 slope of t he power 
spect rum. We used the treePM-SPH code GADGET2 jSpringell 
I2OO5I) to simulate the evolution of a cosmological box with side 
length of I/box = 64^^^Mpc. Within this box we identified (in 
a lower-resolution run utilizing 1024'^ particles) the position of a 
model local group tha t closely resembles the real Local Group (cf. 
ILibeskind etalj20I(]|) . This Local Group has then been re-sampled 
with 64 times higher mass resolution in a region oi2h^^ Mpc about 
its centre giving a nominal resolution equivalent to 4096'' particles, 
resulting in a mass resolution of ttidm = 2.5 x 10^/i~^Mq for the 
dark matter and mgas ~ 4.4 x IO'^/i^^Mq for the gas particles. 
The force resolution is 0.15/i^^ kpc. 

For this particular study we focus on the gas dynam- 
ical smoothed particle hydrodynamics (SPH) simulation, in 
which we follow the feed back and star formation rules of 
ISpringel & Hemguisll ( l2003h : the interstellar medium (ISM) is 
modelled as a two phase medium composed of hot ambient gas 
and cold gas clouds in pressure equilibrium. The thermodynamic 
properties of the gas are computed in the presence of a uniform but 
evolving ultra-violet cosmic back ground generated from QSOs and 
AGNs and switched on at z = 6 faaardt & Madaull 19961) . Cooling 
rates are calculated from a mixture of a primordial plasma compo- 
sition. No metal dependent cooling is assumed, although the gas 
is metal enriched due to supernovae explosions. Molecular cooling 
below lO^K is also ignored. Cold gas cloud formation by thermal 
instability, star formation, the evaporation of gas clouds, and the 
heating of ambient gas by supernova driven winds are assumed to 
all occur simultaneously. 

To give a visual impression of our simulated Local Group we 
present in Fig.[T]a colour-coded density map of the large-scale envi- 
ronment as well as a zoom into the Local Group of our constrained 
simulations. For more details we refer the reader to the CLUES 
web site or the aforementioned papers. But please note that the ac- 
tual simulation is not up for debate here: the prime objective of 
this work is to compare objects finders and hence any reasonable 
simulation of galaxy formation might serve this purpose. 
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in lGill et alj | |2004) as well as lKnollmann & Knebel ( l200^ . The ini- 
tial particle lists are obtained by a rather elaborate scheme: for each 
subhalo the distance to its nearest more massive (sub-)halo is cal- 
culated and all particles within a sphere of radius half this distance 
are considered prospective subhalo constituents. This list is then 
pruned by an iterative unbinding procedure using the (fixed) sub- 
halo centre as given by the local density peak determined from an 
adaptive mesh refinement hierarchy. For more details we refer the 
reader to the aforementioned code description papers as well as the 
online documentation. 

Gas Treatment For the data given here AHF version v 1.0-021 
was run in two configurations: one where the thermal energy of 
the gas entered the unbinding procedure via 



Egas = 4>+ -V 



(1) 



where egas is the total specific energy (required to be negative for 
bound particles), (f> the gravitational potential, v the gas particle's 
velocity and u — its thermal energy, and a second one 

where this feature was switched off (i.e. u = 0). 



Figure 1. A graphical representation of the simulation data showing both 
the (constrained) large-scale structure of the simulation as well as the 
(zoom) region containing the three main galaxies M31, the MW, and M33. 
In Section O the properties of those three main galaxies will be compared 
whereas Sections |5]and|6]focus on the satellite populations of M31 and the 
MW. Image courtesy CLUES collaboration. 

2.2 Lighting up (Sub-)Haloes 

The stellar populati on synthesis model STARDUST (see 
iDevriendt et al. and references therein for a detailed 

description) has been used to derive luminosities from the stars 
formed in our simulation. This model computes the spectral energy 
distribution from the far-UV to the radio, for an instantaneous star- 
burst of a given mass, age and metalicity. The stellar contribution 
to the to tal flux is calcu lated assuming a Kennicutt initial mass 
function jKennicuttll 19981) . 



3 THE CODES 

Before (briefly) presenting the particulars of the codes, it is impor- 
tant to state that all of the participating codes are at various stages 
with respect to their capabilities of handling simulation data that in- 
cludes not only (equal) mass dark matter particles, but also star and 
gas particles. For instance, SUBFIND with its careful treatment 
of the gas thermal properties, is the most advanced with respect to 
the simulation data at hand, whereas on the other hand ROCKSTAR 
required substantial post-processing as it treats all particles on a 
equal footing at present. 

3.1 AHF (Knebe) 

The halo finder AHljfl (AMIGA Halo Finder) is a spherical over- 
density finder that identifies (isolated and sub-)haloes as described 

^ AHF is freely available from http : / / www . popia . f t . uam . es / AHF 



3.2 JUMP-D (Casado & Dominguez-Tenreiro) 

JUMP-D is a galaxy finder and not a sub-halo finder and hence 
is treated differently than the other finders in this work. It aims at 
finding and measuring central and satellite galaxies within given 
host haloes, i.e. baryonic substructure objects within a sphere of 
given radius about the centre of the host. To this extent the 
stellar and gas mass profiles are searched for jumps (and hence the 
name) in the three-dimensional cumulative mass profiles from the 
host halo centre out to the limiting radius (i e- usually the host 
halo's virial radius). The jump detection criterion is based on the 
detection of changes in the first and second derivatives of the re- 
spective mass profiles in the r, 9 and 4> variables at the substructure 
locations corresponding to the humps they cause. For the stellar 
object, the jump in the stellar mass profile is used as a first satel- 
lite detection (i.e., location and velocity), that is later on refined by 
searching for maxima in 6-dimensional phase-space within an al- 
lowance region about that first center, returning the object stellar 
sizes Tstar as well. The jumps in the gas profile are then matched 
to the stellar objects and gas particles inside a spherical region de- 
fined by the radial extend of the gas jump (rgas) are then associated 
to the stellar object. Note that for the detection of the jumps only 
cold gas is considered. 

Please note that the approach of JUMP-D is substantially dif- 
ferent to halo finders in general. The code only locates a bary- 
onic object ("galaxy") without considering the dark matter. To 
this extent JUMP-D cannot be subjected to the common post- 
processing pipeline (cf. Section |T6t when it comes to subhaloes 
as that pipeline heavily relies on the embedding of satellite galax- 
ies within dark matter subhaloes. Therefore, JUMP-D only con- 
tributes to those plots where the identities of such dark matter sub- 
halo particles is not required. The situation is different for field/host 
haloes where it is more straight forward to identify the surround- 
ing dark matter halo and its virial radius. In that regards, JUMP-D 
has been passed through the pipeline for the results presented in 
Section|4] 

Gas Treatment For the satellite galaxies we take advantage of the 
bi-phase nature of the gas in this particular simulation, and just 
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consider the cold gas when detecting jumps in the gas profile. To 
determine the cut in internal energy separating both the cold and hot 
gas regimes, the thermal energy distribution of gas particles within 
the host virial radii is studied and its minimum is used to separate 
the two phases and reject the hot gas from the object, respectively. 
However, this is only applied to the satellite galaxies and not the 
three main objects where all gas irrespective of its phase is added. 

3.3 ROCKSTAR (Behroozi) 

ROCKSTAR (Robust Overdensity Calculation using K-Space 
Topologically Adaptive Refinement) is a phase-space halo 
finder designed to m aximise halo consistency across timesteps 
teehroozi et"ai]|201lh . The algorithm first selects particle groups 
with a 3D Friends-of-Friends (FOF) variant with a very large link- 
ing length (6 = 0.28). For each main FOF group, Rockstar builds 
a hierarchy of FOF subgroups in phase-space by progressively and 
adaptively reducing the linking length, so that a tunable fraction 
(70 per cent, for this analysis) of the particles are captured within 
each subgroup as compared to the immediate parent group. When 
this is complete, Rockstar converts FOF subgroups into seed haloes 
beginning at the deepest level of the hierarchy. If a particular group 
has multiple subgroups, then particles are assigned to the subgroup 
seed haloes based on their phase-space proximity. This process is 
repeated at all levels of the hierarchy until all particles in the base 
FOF group have been assigned to haloes. Unbinding is performed 
using the full particle potentials; halo centres and velocities are cal- 
culated in a small region close to the phase-space density maxi- 
mum. Rockstar is a massively parallel code (hybrid OpenMP/MPI 
style). 

Gas Treatment As already mentioned, initially ROCKSTAR pro- 
vided a particle ID list that was based upon an analysis that treated 
each particle equally (without even taking into account the differ- 
ence in mass between the phases). In the common post-processing 
stage (see Section [X6l below for more details) an unbinding pro- 
cedure was then applied outside of ROCKSTAR. However, this un- 
binding (deliberately) did not take into account the thermal energy 
of the gas, but only considered the kinetic and potential energy of 
all the particles. 

3.4 STF (Elahi) 

The STructure F inder is a parallel hybrid OpenMP/MPI code 
( lElahi et al]|201ll , STF) that identifies objects by utilizing the fact 
that dynamically distinct substructures in a halo will have a local 
velocity distribution that differs significantly from the mean, i.e. 
smooth background halo. This method consists of two main steps, 
identifying particles that appear dynamically distinct and linking 
this outlier population using a Friends-of-Friends-like approach. 
Since this approach is capable of not only finding subhaloes, but 
tidal streams surrounding subhaloes as well as tidal streams from 
completely disrupted subhaloes, for this analysis we also ensure 
that a group is self-bound. Particles which are gravitationally un- 
bound to a candidate subhalo are discarded until a fully self -bound 
object is obtained or the object consists of fewer than 20 particles, 
at which point the group is removed. 

Gas Treatment For this study, we treat each particle type sepa- 
rately, first identifying dark matter (sub)haloes. Gas and star par- 
ticles are then associated with the closest dark matter particle in 



phase-space belonging to a (sub)halo. Further, in order to highlight 
differences in how particles can be treated, we do not pass bary- 
onic particles through an unbinding routine. Part of the motivation 
for this comes from observational studies. It is quite difficult to de- 
termine from observations alone whether gas or stellar structures 
are bound to a galaxy. Instead, spatially coincident gas and stellar 
structures lying within the same line-of-sight velocity window are 
generally assumed to be dynamically related. Hence in this spirit, 
we do not attempt to determine whether a gas or star particle is 
bound to a dark matter (sub)halo, nor do we attempt to account for 
the thermal properties of the gas. 

We note that multiple masses do not effect the local velocity 
distribution function estimate though they do affect the background 
velocity distribution estimate, specifically the centre-of-mass of the 
coarse grain cells used to estimate the local velocity distribution, 
are account for when unbinding the dark matter haloes). 

3.5 SUBFIND (Dolag & Springel) 

SUBFIND identifies substructures as locally overdense, gravita- 
tionally bound groups of particles. Starting with a halo identified 
through the Friends-of-Friends algorithm, a local density is esti- 
mated for each particle with adaptive kernel estimation using a pre- 
scribed number of smoothing neighbours. Starting from isolated 
density peaks, additional particles are added in sequence of de- 
creasing density. Whenever a saddle point in the global density 
field is reached that connects two disjoint overdense regions, the 
smaller structure is treated as a substructure candidate, followed by 
merging the two regions. All substructure candidates are subjected 
to an iterative unbinding procedure with a tree-based calculation 
of the potential. The SUBFIND algorithm is discussed in full in 
ISpringeletallfeOOlb and its extension to dis sipative hydrodynam - 
ical simulations that include star formation in iDolagetaLl l laDOgl) . 

Gas Treatment Within SUBFIND the total density at the position 
of each particle is estimated by summing up the contribution of all 
the different particle species. In contrast to applying the SPH for- 
malism to all particles at the same time, this allows a much more 
fine structured density field to be obtained if the particles of the 
different species have significantly different spatial distributions, 
as is usually the case for star particles. We use all the particles ini- 
tially belonging to the substructure candidate in order to evaluate 
the gravitational potential. For gas particles, we also take the in- 
ternal thermal energy into account in the gravitational unbinding 
procedure. If the number of bound particles left is larger than 20 
and the object contains at least one star or one dark matter particle, 
we register the substructure as genuine sub-halo. 

3.6 Common Post-Processing Pipeline 

As was the case for the previous comparison projects we again sub- 
jected all halo finder (excpet for JUMP-D satellite galaxies) results 
to a common post-processing pipeline. Code representatives were 
simply asked to return the particle IDs of those particles they con- 
sidered to belong to a given object (may that be a host or sub-halo). 
To this extent, halo positions are iteratively determined centre-of- 
masses using the innermost 5 per cent of particles (limiting the po- 
sitional shift per iteration to the force resolution of the simulation, 
i.e. O.I5/i~^ kpc), the bulk velocity is the mean velocity of all par- 
ticles, the mass corresponds to A/200 (i-e. the mass enclosed in a 
sphere so that the mean density inside the sphere equals 200 times 
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Table 1. Some general properties of the three main galaxies. Masses are measured in 10^^/i~^Mq, velocities in km/sec, magnitudes are the Johnson V-band 
magnitude, and the baryon fractions ai'e given as fx = Mx /A/tot where fx can be either /g (gas mass fraction) or fs (stellar mass fraction). Note that the 
total baryon fraction is the sum fb = fg + fs- 



code Mass Vmax Johnson hiv gas mass fraction fg stellar mass fraction fs 
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Figure 2. Density profile of the total (first column), dark (second column), gas (third column), and stellar (fourth column) matter for M31 (upper row), the 
MW (middle row), and M33 (bottom row). The vertical line (coinciding with the 'kink') in the gas profile at 8, 6, and \0h~^ kpc, respectively, indicates the 
transition from the central baryon dominated region to the halo and will be used for the bulge-disk decomposition in Section l431 



the critical density of the Universe), and Knax is the peak value 
of the circular velocity curve. This approach entails that any scat- 
ter reported here will be a lower limit as all halo properties are 
defined in an identical manner. The star particles and their metalic- 
ities and ages are used to calculate the luminosities and magnitudes 
(cf. Section l2l2l for the model). Note that we primarily use magni- 
tude tliroughout this work, liberally referring to it as "luminosity" 
in places. 

We also need to mention that the mode of operation of this 
pipeline biases results towards galaxies residing in dark matter 
haloes: both the definition of the object's edge and the fact that 
we require Vmax to be calculated relies on the fact that the bound 



dark matter particles within such a subhalo have been identified. 
Therefore - as mentioned before - the JUMP-D satellite outputs 
are not best suited to be passed through this pipeline. 



4 GALAXY COMPARISON 

The present work would not rightfully be called a "Galaxy Com- 
parison Project" if we were not to study the properties of the (bary- 
onic component of the) dominant galaxies in our constrained Local 
Group simulation, i.e. the MW, M31, and M33. We therefore begin 
by comparing the particulars of these objects. 
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Figure 3. Same as Fig.|2]but this time showing the relative difference to the mean of all density profiles. 



4.1 General Properties 

In Table [T] we provide a summary of the most fundamental prop- 
erties of the three main galaxies in our simulation, i.e. we list the 
mass, Vmax, luminosity (as characterised by the Johnson V-band 
magnitude), and baryon fractions fx = Mx/Mtot where fx can 
be either fg (gas mass fraction) or fa (stellar mass fraction). Note 
that the total baryon fraction is the sum fb = fg + fs - There is a 
rather excellent agreement across the halo finders, with the scatter 
in basically any quantity being smaller than 10 per cent, which was 
the deviation found in our previous dark matter only comparisons. 
The remaining differences seen in the mass are readily explained 
by the fact that some codes (e.g. AHF) consider subhaloes to con- 
tribute to the mass of their respective host while others (e.g. STF 
and SUBFIND) exclude the subhalo masses. 



4.2 Mass Distributions 

To explore any differences in more detail we now decompose the 
mass profiles of the three galaxies. The result is presented in Fig.|2] 
where we show the total matter density profile (left column) and 
its breakdown into the various components, i.e. dark matter (2nd 
column), gas (3rd column), and stars (rightmost column). The rows 
represent the MW (upper row), M31 (middle row), and M33 (bot- 
tom row). 

All finders agree outstandingly well for all three components 
(and hence the total matter, too). Despite their different treatment 
of the gas thermal energy, the resulting mass distributions through- 
out the three galaxies remain essentially indistinguishable. When 
comparing the profiles to the numbers in Table [T] one may wonder 
about M31 and why it is that while the masses differ by 30 per 
cent the density profiles agree so remarkably well. This is read- 
ily explained by differences in the virial radii of approximately 10 
per cent, something not clearly visible in Fig. [2] due to the loga- 
rithmic X-axis. In that regard, note that all the curves in the fig- 
ure terminate at the respective virial radius and go slighter closer 
to the centre than the nominal force resolution of the simulation, 
i.e. 0.15/i~^ kpc. These variations in mass and radius are partly 
due to different treatment of subhalo masses: some codes include 
substructure in the host mass (e.g. AHF) whereas other do not 
(e.g. SUBFIND). And this is particularly important for the sim- 



ulated M31 as it contains one large group of satellites contribut- 
ing approx. 15 per cent of the total host mass (cf. Sec. 4.2.1 in 
iKlimentowski et al]l201(l where this particular object is discussed 
in detail). However, the mass and radii differences are also affected 
by disparate particle collection and unbinding procedures (Knebe et 
al., in preparation) accounting for the remaining variations across 
codes. 

Another interesting aspect of Fig.[2]is the kink in the gas den- 
sity profile apparent at approximately 6-10/i^^ kpc, corresponding 
to roughly 10 per cent of the virial radius (with no difference across 
finders again). This peak represents the transition from the central 
baryon dominated region to the halo and will be used in the follow- 
ing Sub-Section for the bulge-disk decomposition. 

To better gauge the differences across finders we addition- 
ally plot in Fig. [3] the fractional differences between the profile of 
each finder and the respective mean profile when averaging over 
all codes. This figure shows that the residual scatter of the density 
profiles is at most just 10 per cent about the mean. 



4.3 Dynamical Bulge-Disk Decomposition 

A dynamical bulge-disk decomposition for all three central galax- 
ies in our simulation has been performed. We decompose the galax- 
ies using two different techniques, one for the gas and one for the 
stars. 

For th e star particles we apply the metho d proposed by 
lAbadi et all (l2003h. used bvlOkamoto et alj | |2005|) , and improved 



upon in Domenech-Moral et alj ( l20I2h . The method is to essen- 
tially compare the ^-component of the angular momentum (J^) of 
each particle with the angular momentum of a circular orbit with 
the same energy, Jc(E ). The first step is to thus define a ^-direction 
which, as suggested bv lAbadi et alj ( l2003l) may be chosen to be the 
total (specific) angular momentum of all the star particles inside the 
"luminous radius" of the galaxy (defined here as within 6, 8, and 
lOh^^ kpc for the MW, M31, and M33, respectively; see the kink 
in the gas profiles visible in Fig. |2j. The choice of a z-direction 
in this manner ensures that the majority of the particles will have 
positive values of . The potential and kinetic energy for each star 
and gas particle within this luminous radius of the galaxy is then 
calculated. Since circular orbits maximise angular momentum, we 
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Figure 4. A dynamical bulge - disc decomposition for stars (top) and gas (bottom) particles in the MW (left), M31 (middle) and M33 (right). We compare each 
star particle's angular momentum in the z-direction to the angular momentum of a circular orbit of the same energy and show on the j/-axis the probability 
density distribution P(e), for a given e = Jz/Jc{E). For gas particles we compare the angular momentum in the z-direction with the angular momentum of 
a ckcular orbit at that radius and show on the y-axis the probability density distribution P(e), for a given e = Jz/ Jcirc- 



assume Jc{E) for a particle of energy E is simply the maximum 
value of Jz of all particles with that energy. Note that the limits of 
e = Jz/Jc{E), i.e. -1 and 1, correspond to counter- and co-rotating 
particles on circular orbits. Spheroidal systems with very little net 
rotation are thus represented by a Gaussian centred on 0. A fully 
rotationally supported thin disk on the other hand is represented by 
values of Jz/Jc{E) ~ 1. Thick disks not fully supported by their 
rotation can be seen as distributions cen tred on positive values such 
as Jz/Jc{E) ^ 0.5 l lAbadi et al.ll2003l) . 

The resulting distributions can be viewed in the upper panels 
of Fig. |4] All three galaxies contain both co- and counter rotating 
populations. The MW comprises a bulge and a (co-rotating) thick 
disk that peaks at e « 0.8. No thin disc is seen here, although ow- 
ing to their relatively high values of e, many particles in the thick 
disc are kinematically cold. M31 on the other hand appears as a 
featureless spheroid whose net angular momentum can be seen as 
a fat tail towards positive values of e. M33 is a more model galaxy 
- it contains a well defined rotationally supported thin disc as well 
as a central bulge and a small, kinematically hotter, thick disc com- 
ponent. 

While all this is certainly interesting in terms of galaxy forma- 
tion, the most important part to note - with respect to this compar- 
ison project - is that all halo finders behave similarly when exam- 
ining this stellar decomposition. 

For the decomposition of the gas particles we use a slightly 
different method because some of the halo finders include the in- 



ternal energy in their unbinding procedure, while others do not. Un- 
like the stars, it is thus unfair to decompose the gas distribution on 
energetic grounds. Consider a co-rotating stellar and gaseous disc. 
Some halo finders will use the internal energy of gas particles to un- 
bind them from this coherent motion and thus exclude them from 
the dynamical decomposition described above. In this case, while 
all halo finders will easily pick out star particles that are on circular 
orbits, not all gas particles on circular orbits will be included. The 
results of such a dynamical decomposition comparison are thus un- 
reliable. 

In o rder to overcome this inco nsistency we use a method pro- 
posed bv lScannapieco et"al] 1 I2OIOI) . As with the stars, a z-direction 
is defined to be the total angular momentum of all gas particles 
within 6, 8, and lOh^^ kpc for the MW, M31, and M33, respec- 
tively. We then calculate the ratio of e = Jz/Jdic where Jcirc is 
the angular momentum of a circular orbit at the same radius, i.e. 

Jcirc — r X ij'^^^^^. In this way we bypass the need to use 
the potential and kinetic energies since, as mentioned earlier, gas 
particles also have thermal energies. When this ratio is equal to 
unity, a gas particle is moving on a circular orbit (with respect to 
the 2-direction). Note that this ratio is in principle unbounded from 
(—00,00). In practice however, few (~ 1 per cent gas particles 
counter rotate and thus a lower limit of zero is applicable. 

The decomposition of the gas particles is shown in the bot- 
tom panels of Fig. |4] All gas distributions exhibit peaks at e = 1, 
indicating the presence of gaseous discs. M31 has a wider distri- 



© 2010 RAS, MNRAS OOO.rTlfTsl 



8 Knebe et. al 



Table 2. Number of subhaloes (compliant with our selection criterion, see 
text for details) found by each participating halo finder. For JUMP-D we 
list the number of baryonic objects found. 
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bution implying a thicker, dynamically hotter component (as also 
suggested by the stellar decomposition above it). It is interesting to 
note that in M31 AHF and RoCKSTAR identify a small local min- 
imum around e = 1. This dip (with different maxima to the left 
and the right) is interpreted as a subtle warping of the gaseous disc. 
The astute observer will also notice a small but non-negligible tail 
towards zero in the MW and M31 distributions; these are gas par- 
ticles whose angular momentum in the z-direction is smaller than 
that of a circular orbit of this radius and who are thus on either 
radial or inclined orbits. 

The minor discrepancies between halo finders for the gas com- 
ponent are certainly due to the differences in the handling of gas 
particles: for instance, AHF as well as SUBFIND read the gas 
thermal energy taking it into account during the unbinding proce- 
dure; all other codes do not consider the thermal energy at that stage 
(STF and RoCKSTAR) or do not feature an unbinding procedure 
(JUMP-D). 



5 SATELLITE COMPARISON 

We now turn to taking a closer look at the satellite populations of 
the two prominent host haloes, i.e. all subsequent plots are based 
upon a combined sample of MW and M3 1 subhaloes that reside in- 
side a sphere of radius 300 kpc centred on their respective host and 
that contain at least one star particle. Note that this radius is slightly 
larger that the virial radii of the two hosts ( ~ 200 kpc) and moti - 
vated by the observationally inferred value l lWatkins et allbOlOh . 
We further impose a lower limit for V^nax lOkm/sec as this cor- 
responds to the completeness level (cf. Fig.[5]below). We list in Ta- 
ble |2] the number of subhaloes compliant with these criteria found 
by each of the participating finders.; and for JUMP-D we give the 
number of baryonic objects found. Further note that JUMP-D is in 
the special situation of not considering subhalo finding and hence is 
not suited to be subjected to the common post-processing pipeline. 
Therefore, JUMP-D only contributes to the plots that merely re- 
quire knowledge of the baryonic component. This data has been 
directly provided by the finder itself. 



5.1 Vmax Function 

We start with a quantity that has caught a lot of atten- 
tion recently, especially for the Local Group/Milky Way dwarf 

the peak of the circular velocity 



curve Krav (e.g.|Bovlan-Kolchin et alj20Ilb 


;ldiCintio etalj201ll; 


Bovlan-Kolchin et al.ll201 lal; Vera-Ciro et al. 


2OI2I; Di Cintio et all 


20I2|). It has. for instance, been shown bv Idi Cintio et alj J201 ll) 



that the inclusion of baryonic physics into the simulation can lead to 
variations in Vmax (and -Rmax) when being compared to dark matter 
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Figure 5. Cumulative Vmax function for haloes closer than 300 kpc to ei- 
ther of the two simulated hosts M3 1 and MW that at least contain one star 
particle. 
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Figure 6. The relation between the peak of the circular velocity curve Vmax 
and its position -Rmax for subhaloes inside a 300 kpc sphere about either 
M3 1 or the MW. The thin solid lines dehneat e the la confidence interval of 
the observed bright M W dSph galaxies (as in iBovlan-Kolchin et alj20TTbl : 
Idi Cintio et alj201 ih . 



only simulations; in particular, (sub-)haloes with high baryon frac- 
tions tend to experience adiabatic contraction while objects with 
lower baryon content move towards l ower- Kn ax/higher- -R max val- 
ues, possibly due to mass outflows jNavarro et al.lll99^ or ran- 
dom bulk motion of gas "heating" the central matter distribution 
jMashchenko et alj2006l) . 

In Fig. |5] we compare the cumulative distribution of Vinax 
found by all our finders. We find excellent agreement. Please note 
that this is the only plot that has not been subjected to the lower 
Vmax hmit. 



5.2 -Rn 



■Vm 



relation 



As already outlined above when motivating the consideration of 
the 'Knax function in the first place, there has been some substan- 
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tial debate about the failure of the most massive subhaloes found 
in simulations of Milky Way type d ark matter host haloes to host 
the most luminous observed dSph's iBovl an-Kolch in et al. 201 Ibl; 



number of star particles 
10 100 1000 10000 



diCintio etai]|201ll: 



2012l : lDi Cintio et al 



Bovlan-Kolchin et alji201 la> . iVera-Ciro et alj 
2012h . To verify whether different halo find- 



ers will contribute differently to this discussion we present in Fig.|6] 
the usual plot confronting i?max with Vmax- The two solid lines 
delimit the la confidence interval of the observed bright Milky 
Way dwarf spheroidal g a laxies , as defined in lBovlan-Kolchin et alJ 
( l2011bl) : ldi Cintio et alj ( l2oTll . please refer to these references for 
more details). 

We again find remarkable agreement across the participating 
codes, even though there appears to be some marginal scatter across 
Rmax values. This can be attributed to the fact that the determina- 
tion of Knax is certainly more stable for cosmological dark matter 
haloes due to the flatness of the circular velocity curve; and this 
flatness is obviously also responsible for (slight) variations in the 
position of the peak velocity. In that regard, it should be noted that 
in order to best calculate the position of the peak, i.e. 7?max, in the 
circular velocity curve we decided to only use the dark matter for 
this purpose, i.e. i?max is defined as the position of the maximum 
of i\fDM(< r)/r\ however, V^ax = GMtot(< -Rmax)/-Rmax ob- 
viously takes into account all matter. 



5.3 Baryonic Mass 

In Fig.|7]we show the cumulative stellar (top) and gas mass (bot- 
tom) function of all our subhaloes. The number of objects in the 
stellar mass plot corresponds to the one given in Table |2] But we 
also acknowledge that this number is smaller than the total number 
of "subhaloes with stars" found by the majority of finders as can be 
seen in Fig.|5]due to the now applied Kiiax cut. 

While the stellar masses of the objects agree remarkably well 
(note that JUMP-D returns the stellar masses of the objects within 
rstar, and hence their slightly smaller values in some cases), we no- 
tice quite substantial differences in the gas masses. While we return 
to the differences found for the gas fractions in subhaloes in more 
detail later, we bear in mind that the actual number of subhaloes 
containing any gas is rather low for the majority of finders (note 
that JUMP-D returns the cold gas masses within rgas only when 
cold gas exist at the satellite center). STF, which does not perform 
an unbinding procedure for the gas, finds a substantial number of 
satellites containing gas particles down to the limit of one gas parti- 
cle. The mass of a gas particle is rrigas = 4.4 x lO^/i^^M© which 
corresponds to the lower limit of the x -axis. 

While the agreement of the (stellar) masses appears to indi- 
cate that the objects found by the suite of finders used here are 
the same, this is not necessarily the case given their different modi 
of operations. To shed some more light o n the cross-correlati on of 
objects (in a statistical sense) we follow lOnions et aL I ( l2012h and 
plot in Fig. [8] the cumulative stellar mass in objects as a func- 
tion of distance from the centre of the host halo (normalized to 
the host radius). We stack the data from the MW and M3I. The 
agreement proves that the satellite galaxies not only have similar 
masses but also live at approximately the same distance from the 
central galaxy. There are, however, a few (massive) objects found 
by JUMP-D and STF close to the centre not picked up by the other 
finders. There is another interesting feature in this plot, namely that 
JUMP-D is slightly below the other finders at larger radius, which 
is possibly due to ignoring the underlying dark matter, so that ei- 
ther the individual galaxies have less stars or some with a low stellar 
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Figure 7. Cumulative stellar (top) and gas (bottom) mass functions for the 
same objects already shown in Fig.|6] The upper x-axis of each panel pro- 
vides a translation of mass into number of particles. 
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Figure 8. Cumulative stellar mass in objects as a function of distance to the 
(noraialized) host centre. 



© 2010 RAS, MNRAS OQO.rTlfTSl 



10 Knebe et. al 




10 100 -5 -10 -15 -20 

I'm., [km/s] [mag] 



Figure 9. Baryon fraction /[, vs. Vm 



content are missed because without the dark matter particles they 
are too small to be assigned as reliable substructures. 

So far we have separated stellar and gas mass when studying 
baryons. But we now briefly turn to the baryon fraction, investigat- 
ing its relation with Vmax in Fig.|9] Again, all codes show the same 
tendancy for more massive subhaloes to also contain more baryons 
and it is difficult to decipher any differences (despite our knowl- 
edge that there are marked differences in the gas content), unless 
we move to a comparison of objects on a one-to-one basis to be un- 
dertaken in Section|6] However, we note that STF does not find as 
many low-/;, objects (for Knax < 20 km/sec) as the other finders 
due to the fact that it keeps more gas in its subhaloes as seen in (the 
bottom panel of) Fig.|7] Bear in mind again that STF did not pass 
the gas particles through an unbinding routine, hence the expected 
large gas fractions at all mass scales. 



5.4 Luminosity Function 

In Fig. [To] we now present the Johnson V-Band lu minosity func- 
tion as well as the observa ti onal d ata as taken from lKoposov et al] 
( l2008l) and iMaccio et al. respectively (thin solid line, 

referred to as "Maccio sample"): these data are a combina- 
tion of the volume c orrected MW satellite luminosity function 
Koposov et al. 20081) augment ed with information from iMated 



fvoposov et al.l UUUel) augment ed witn iniormation trom IMated 
19981) and Maccio et al.l ( I2OI0I) kindly provided to us by Andrea 



Maccio (personal communication). And even though all luminos- 
ity functions agree with the Maccio sample rather well, we stress 
that we included the observational data merely as a reference to 
guide the eye. It is not our prime objective to reproduce the MW 
and/or M3 1 luminosity function of satellite galaxies with our simu- 
lation. We find again excellent agreement between the finders. But 
his comes at no surprise as we found negligible variations in the 
stellar mass of subhaloes in (the upper panel of) Fig. |7] already. 



6 SATELLITE CROSS-COMPARISON 

While the general agreement for the satellite populations between 
codes is rather outstanding - at least when it come to the fun- 
damental characteristics such as Vmax and luminosity My - we 
now go one step further and directly compare properties of objects 



Figure 10. The luminosity function of subhaloes in the Johnson V-Band. 
The "Maccio" observational data (thin solid line) is a combinati on of 
the volume coiTected MW luminosity function Koposov et al. ( 2008j) aug- 
mented with information from Mateo ( 1998) and Maccio et al. |2010[) under 
the assumption of an NFW-like radial distributions of satellites. Note that 
the comparison to the observational data is not the prime target of this study 
and only serves as a reference. 



on an individual basis. To this extent we use the same matching 
criterion (based upon the uniqueness of particle IDs) successfully 
used to construct merger trees or cross-cor r elate different simula- 
tions before (e.g 



iKnebe et aLllioTl 



Klimentowski et aljboid : iLibeskind et"ai]|20ld: 



y) by facilitating a tool that comes with the AHF 
package called MergerTree. Even though originally designed to 
identify corresponding objects in the same simulation at different 
redshifts it can equally be applied to find "sister objects" in an anal- 
ysis of the same snapshot done with a different halo finder. 

Note that distribution functions (as well as scatter plots) previ- 
ously presented will not only suffer from differences in individual 
halo properties but also encode the fact that some finders may have 
identified different numbers of objects (cf. Table|2j. To circumvent 
this we aim at directly comparing quantities on a halo-to-halo basis 
and move from general distribution functions and their variations 
as discussed above to cross-comparing satellites. The plots in the 
following Sub-Sections |6.2| through to |6.4| now all follow the same 
logic: the x-axis shows the median of the haloes' Vmax, whereas the 
y-axis gives the normalised difference between the lower and upper 
percentiles equivalent to the 2nd and 3rd ranked of the distribution 
across all halo finders (JUMP-D is not included in these plots as 
no dark matter halo required for the calculation of Vmax has been 
provided). We deliberately chose to use medians and percentiles as 
the distribution of properties across finders might be non-Gaussian 
and at times biased by just one outlier. Further note that not all sub- 
haloes have an error greater than zero: those for which the differ- 
ence in the 3rd and 2nd ranked finder is zero have been artificially 
placed at the bottom of each panel at a value of 1.2 x 10^'* yet 
entered with their correct value of zero into the calculation of the 
error percentages given below in the respective Sub-Section. We 
further indicate the 1 per cent and 10 per cent error by a horizontal 
dashed line in each plot and provide the fraction of those objects 
below that error as a legend. 
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6.1 A Common Set 

For a given subhalo in one catalogue we locate the subhalo which 
contains both the greatest fraction of its particles and which is clos- 
est to it in mass. Specifically, for each sister candidate i we calcu- 
late a merit rrii, defined as 



rrii 



2 

^shared,] 



^reference ^sister, z 



(2) 



where Tisharod.i is the number of particles shared between both sub- 
haloes, riaister,! is the total number of particles in the sister candi- 
date and nrefcrcncc IS the total number of particles in the reference 
subhalo. We then identify the sister as that halo with the largest 
value of rrii. Typically, successful matches result in rrn ^ 0.5 and 
we assume to have failed to find a sister if rrii < 0.2 for a given 
subhalo[3 

By restricting ourselves to the set of objects found by every 
halo finder we are able to directly compare the properties of the 
same object across all finders. Even though AHF found the most 
haloes compliant with our initial selection criterion and we chose to 
use its catalogue as the basis, we confirm that every subhalo found 
by AHF has an actual counterpart in all other analyses. However, 
some of these counterparts do not fully lie within the required 300 
kpc-sphere and hence were rejected from the previous analysis. Us- 
ing them now will increase our statistics without obscuring the re- 
sults: the aim is to compare the same objects across different finders 
and they certainly serve this purpose. 



6.2 Knax&i? max 

While we have seen that the general relation between the codes for 
subhaloes compliant with our selection criterion is rather remark- 
able. Fig. [TT] further confirms that also the one-to-one scatter of 
cross-identified objects is small: inspecting the differences in the 
peak of the circular velocity curve Vmax and its location -Rmax, we 
find the overall median values of the error to be < 1 per cent and 
2 per cent for Vmax and i?max, respectively. It is evident that Knax 
gives better agreement as it is a more stable property than i?max 
due to the flatness of the circular velocity which is responsible for 
at least a part of the scatter in -Rmax- 



6.3 Baryons 

The comparison of the gas mass function in Fig.|7]already revealed 
differences while the stellar mass function showed excellent agree- 
ment. It therefore comes as no surprise that we find a median error 
of 86 per cent for the gas mass Mgas (considering only subhaloes 
that actually contain gas) and < 1 per cent for the stellar mass A/* 
in Fig. [12] As already noted, some subhaloes may be missing from 
the plot if the difference between their 3rd and 2nd ranked find- 
ers is zero. This is in fact the case for some of the objects entering 
the stellar and gas mass comparison: 47 (out of the 75 success- 
fully cross-matched subhaloes) show identical stellar masses and 
63 have no gas particles at all. In that regard, the actual median 
error of the gas mass is per cent (while the stellar mass value 
remains unaffected), but we decided to rather report the deviation 

We found these criteria to give the most reliable merger trees for sub- 



haloes (e.g. Klimentowski et al. 2010; Libeskind et al. 2011; Knebe et al.' 
[201 1 ; di Cintio et al. 2011 ; Libeskind et al. 2011 ; Di Cinti o et aU2012i) and 
hence decided to apply them also for the cross-correlation. 
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Figure 11. Relative error in Vmax (top) and -Rmax (bottom). The horizontal 
dashed line simply indicates an error of 1 per cent for reference. Subhaloes 
appealing at an error value of 1.2 X 10^^ are the ones that actually have 
zero values. 



only considering subhaloes containing gas and seen in the plot, re- 
spectively. 

How does the error in baryons with the simultaneous agree- 
ment in stellar mass translate into differences in the baryon frac- 
tion of subhaloes? This can be verified in Fig.[T3]where we find an 
overall scatter of 5 per cent. This leads to the conclusion that the 
stellar mass content is by far the most dominant and important for 
the baryon fraction and the substantial differences in gas mass only 
affect its value marginally. 

6.4 Luminosity 

While the agreement between codes for the V-band magnitude was 
rather excellent - as verified by the luminosity function presented 
in Fig. [To]- will we find similar consistency when comparing ob- 
jects on an individual basis? This can be verified in Fig. 1 141 where 
we find a general median error of < 1 per cent (noting again that 
47 out of 75 cross-matched subhaloes contain an equal amount of 
star particles). One thing to note here is that the errors for Mv are 
a fair bit smaller than the differences found for the stellar masses, 
even though luminosities are directly related to the star particles 
(and their properties like age and metallically). But while the trans- 
formation from the upper panel of Fig. [12] to Fig. [14] is most cer- 
tainly non-linear, the general features (and "satellite placement" in 
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Figure 12. Same as Fig. 1111 but for stellar mass M, (top) and gas mass 
Mgas (bottom). 



the plot) are nevertheless preserved; the objects are only shifted to- 
wards smaller errors entailing that luminosity as measured here via 
Mv is not as sensitive to the stellar content as the mass itself. 



6.5 Mass Distributions 

To further explore the differences - in particular for the gas content 
of subhaloes - we look at the mass distribution of the various com- 
ponents (i.e. dark matter, gas, and stars, as already done for Fig.[2]l 
in the subhaloes found by AHF and cross-matched to all other halo 
finders. In Fig.[T5]we showcase two examples, one massive subhalo 
(top row, containing roughly 20000 particles) and one less massive 
object (bottom row, containing 3500 particles); note that using ob- 
jects with even lower numbers of particles will no longer provide 
credible profiles. From left to right the columns show the total mat- 
ter profile, the dark matter, the gas, and the stellar mass profile (all 
in units of the critical density at redshift z = G). While we can 
again confirm rather excellent agreement for the dark matter and 
stars, we clearly see differences in the gas - as expected from the 
previous plots. In particular, the lower mass subhalo does not fea- 
ture a gas profile at all for the majority of the applied halo finders. 
Note again that all curves terminate at the edge of the object. 



Figure 13. Same as Fig. 1 1 1 1 but for the baryon fraction = (M* 

Mgas)/Mtot. 
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Figure 14. Same as Fig. llll but for the Johnson V-band magnitude My. 

6.6 Gas Treatment - A Test Case 

To better understand the differences seen in the gas properties in 
Figs. [Tineas well as [15] we went back to the simulation data and 
extracted all the particles in a spherical region about the centre of 
the subhalo featured in the lower panel of Fig. [15] We can clearly 
see that the region about the object's centre contains a substantial 
number of gas particles (shown in the upper left panel). But all 
codes featuring a treatment of the gas thermal energy either dur- 
ing or prior to the unbinding (i.e. AHF and SUBFIND) remove 
essentially all gas from the subhalo; the other two finders that do 
not include the thermal energy during the unbinding are left with a 
residual amount of gas. 

When using AHF in the case where the gas thermal energy has 
been ignored, AHF basically considers all gas particles seen in the 
left panel to be part of the subhalo. In contrast, due to their phase- 
/velocity-space nature, both ROCKSTAR and STF consider the ma- 
jority of the gas particles to belong to the background host and keep 
only a small amount of them. For the object considered here the ef- 
fective thermal velocity of each gas particle is always larger than its 
kinetic velocity (not shown here) and hence the grouping in phase- 
or velocity-space will naturally remove (hot) gas whenever a gas 
particle is considered not belonging to it based upon kinetic veloc- 
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Figure 15. Density profile of the total (first column), dark (second column), gas (third column), and stellar (fourth column) matter for a massive (upper row, 
called 'satl47') and a low-mass (lower row, called 'sat623') subhalo. 



ity only. Or put differently, the gas component forming part of the 
background halo is prone to be removed by such finders as they in- 
herently use velocity information when grouping and collecting the 
initial set of particles, whereas configuration space finders only deal 
with velocities (either kinetic or thermal) in a (post-processing) un- 
binding procedure. On a side note, a visual inspection of a larger 
region about this particular sample satellite galaxy indicates that it 
has passed extremely close to its host already and been subjected to 
severe tidal forces; this might also explain why ROCKSTAR appears 
to only have found one side of the object. 

To verify that in actuality all gas particles should have been re- 
moved we finally show in Fig.[T7]the same region again (left panel, 
slightly different projection) alongside the pressure and density of 
the gas particles along the a;-axis (right panel). We can clearly see 
that the gas particles inside the 5h~^ kpc sphere about the centre 
(marked in green) are forming part of the overall background gas 
particles (marked in blue). The stars (red) are only shown as a ref- 
erence. 

In summary, while configuration space finders such as AHF 
and SUBFIND require a proper treatment of the gas' thermal en- 
ergy to remove (unbound) gas particles belonging to the back- 
ground, phase-/velocity-space based finders such as RoCKSTAR 
and STF encode part of this removal into their methodology. How- 
ever, they keep a residual amount of gas. 



7 SUMMARY & CONCLUSIONS 

We set out to compare the performance of object finders for cos- 
mological simulations in the case where the data not only contains 
dark matter but also baryons, i.e. gas and stars. Both types of mat- 
ter obviously cluster more strongly, gas being able to dissipate its 
energy, and star particles forming at the bottom of deep potential 
wells out of the cold gas. While star particles only feel gravity (like 
the dark matter particles), they nevertheless originated from differ- 
ent sites and more compact regions, respectively. Also gas particles 
carry a thermal energy acting like a pressure (given an appropriate 
equation of state). At present, halo finders deal with these subtleties 



differently, with some properly accounting for the gas pressure and 
others ignoring it. Therefore, the fundamental question is, what will 
different halo finders find when it comes to simulations including 
baryonic physics? 

We found remarkable and excellent agreement between dif- 
ferent halo finders when applied to a constrained simulation of the 
Local Group that includes all the relevant baryonic physics; at least 
as far as dark matter and stellar content are concerned. However, 
there appear to be strong discrepancies when it comes to the gas 
and hence baryon fractions. This result may not be too surprising as 
both dark matter and stars are particles primarily only feeling grav- 
ity (even though stars obviously have a different formation history) 
whereas gas is also subjected to hydrodynamical forces. It appears 
that - to properly identify and subsequently analyse and model the 
gas component in cosmological simulations - halo finders need to 
take into account the pressure of the gas in the unbinding procedure, 
especially when dealing with subhaloes containing of order 10,000 
or less particles (in total). But we have also seen that some correct 
accounting of the gas happens in phase-Zvelocity-space based find- 
ers such as RoCKSTAR and STF used here: not all the gas particles 
inside the radius of the subhalo are considered part of the object; 
only a residual number of gas particles are kept. In these finders the 
(hot) gas component is prone to be removed as they inherently use 
velocity information when grouping and collecting the initial set 
of particles, whereas configuration space finders only deal with ve- 
locities (either kinetic or thermal) in a (post-processing) unbinding 
procedure. 

In summary, the proper treatment of the gas certainly also de- 
pends on the actual scientific question to be addressed with the nu- 
merical data; any algorithm used for collecting gas particles should 
carefully consider what type of study they are interested in. For 
instance, for studies where the X-ray emission from (hot) gas in 
galaxy clusters is sought, one clearly would not want to remove the 
majority of the gas (during the unbinding or even before). Further, 
the differences will certainly become more important for simula- 
tions (of subhaloes) with much higher resolution where we expect 
subhaloes to hang on to their own gas; we have seen here that for 
this particular data only a few subhaloes still have gas (cf. Fig.|7]l. 
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Figure 16. Visualisation of a sublialo showing all particles inside a spherical 
region about the identified centre (upper left panel) vs. the actually identi- 
fied particles of the individual halo finder showing all types of particles: gas 
(blue), stars (red), and dark matter (black). Note that the object coincides 
with the lower panel object of Fig. 1151 
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Figure 17. Gas and star particles for the same object as in Fig. 1161 The left 
column shows all gas (blue) and star (red) particles inside a cubical region 
of side length lO/i^^ kpc as well as the gas particles inside a 5h~^ kpc 
sphere (green). The right column shows the pressure (upper) and density 
(lower) for the gas particles along the x-axis using the same colouring 
scheme. 



But consider gas rich objects at high redshift: for these the baryon 
fraction would be dominated by the treatment of the gas and is ex- 
pected to be much larger as stars have not formed yet. This discus- 
sion about the proper treatment of the gas (during the unbinding 
procedure) is akin to debates about unbinding in general: some in- 
vestigations heavily rely on it (e.g. subhalo dynamics) whereas oth- 
ers can live without it (e.g. gravitational lensing). In any case, the 
differences for the three main galaxies investigated here were reas- 
suringly negligible and we conclude that despite subtleties there is 
in general a fair agreement of "galaxy finders". 
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